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Resumen 


Chávez, C., Mota, C., Fuentes, C., € Quevedo, A. 
(enero-febrero, 2018). Modelación bidimensional de la 
infiltración del agua en surcos aplicando el gradiente 
conjugado. Tecnología y Ciencias del Agua, 9(1), 89-100, DOL: 
10.24850 /j-tyca-2018-01-06. 


El conocimiento de la distribución del bulbo de humedad 
durante un evento de riego ayuda a comprender de mejor 
manera la profundidad de mojado. El flujo del agua en suelo 
es modelado con la ecuación de Richards, sin embargo, 
la mayoría de las soluciones requiere de mucho tiempo 
computacional o la solución está acotada a un dominio en 
específico. En este trabajo se resuelve la ecuación de Richards 
2D usando dos esquemas: un esquema implícito en todo el 
perfil y un esquema mixto: explícito en la zona no saturada e 
implícito en la saturada. El esquema implícito se resolvió con 
el gradiente conjugado para problemas no lineales usando 
en la superficie una función sinusoidal. La comparación se 
realizó entre los dos esquemas, y los resultados de perfiles de 
humedecimiento y lámina infiltrada fueron similares, dado 
un criterio de error. El proceso de infiltración se simuló en 
un evento de riego por surcos mediante tirante constante 
en los tres surcos, y en un segundo evento se les impuso 
un tirante diferente a cada uno. Los perfiles de distribución 
de humedad que se obtienen durante el tiempo simulado 
pueden ayudar a dar recomendaciones de las desventajas 
de establecer un tirante diferente en los surcos durante un 
evento de riego. 


Palabras clave: Richards, bidimensional, diferencias finitas, 
riego por surcos. 


Abstract 


Chávez, C., Mota, C., Fuentes, C., € Quevedo, A. (January- 
February, 2018). Modeling infiltration 
of furrow  applying the  conjugate  gradient. Water 
Technology and Sciences (in Spanish), 9(1), 89-100, DOI: 
10.24850/j-tyca-2018-01-06. 


two-dimensional 


Knowledge of the bulb moisture distribution during an irrigation 
event helps to better understand the depth of wet soil. The flow of 
water in soil is modeled with Richards” equation, however, most 
solutions require a lot of computational time or the solution is 
bounded to a specific domain. In this study, the two-dimensional 
Richards' equation is solved using two schemes: one implicit in the 
entire depth and a mixed scheme: explicit in the unsaturated zone and 
implicit in the saturated zone. The implicit scheme was solved with 
the conjugate gradient for non-linear problems using a sinusoidal 
function on the surface. The comparison was performed between the 
two schemes and the results of wetting profiles and infiltrated depth 
were similar, given an error criterion. The infiltration process was 
simulated in a furrow irrigation event through a constant depth in 
the three furrows, and in a second event they were imposed a different 
depth to each one. The moisture distribution profiles obtained during 
the simulated time can help us to give recommendations for the 
disadvantages ofestablishing a different depth in the furrows during 
an irrigation event. 


Keywords: Richards, two dimensional, finite differences, furrow 
irrigation. 
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Introducción 


Las dos modalidades más comunes para el rie- 
go por gravedad son riego por melgas y riego 
por surcos. Estos procesos ocurren, como es 
natural, en el espacio 3D. Cuando la geometría 
del espacio y las condiciones de frontera son 
invariantes respecto a algún eje coordenado, 
el proceso de infiltración puede simplificarse y 
considerarse como un proceso unidimensional; 
este es el caso del riego por melgas, en el cual 
la geometría del espacio es parecida a un para- 
lelepípedo. En el riego por surcos, la geometría 
del proceso es más compleja y requiere que se 
consideren al menos dos variables espaciales 
para su simulación correcta (López-Avendaño, 
Palacios-Vélez, Fuentes-Ruíz, Rendón-Pimentel, 
éz García-Villanueva, 1997). 

La ecuación que representa el movimiento 
del agua en un suelo parcialmente saturado es 
obtenida mediante la combinación del principio 
de continuidad y la ley de Darcy. Si se introduce 
la capacidad específica como la derivada de la 
curva de retención (C(W) = 00/09, la ecuación 
resultante queda expresada en términos de 
una variable 0(W), denominada ecuación de 
Richards (1931), que en su forma bidimensional 
adquiere la forma: 


cy dl. 0 


ot ox AN : AS id 


Ox | 0z OZ OZ 


donde 8 es el contenido de humedad; Y, la 
presión; x, una variable espacial (horizontal); z, 
la profundidad; K, la conductividad hidráulica 
y t es el tiempo. 

Para resolver esta ecuación se incorporan 
las características hidrodinámicas del suelo: la 
curva de retención de humedad y la curva de 
conductividad hidráulica. Las funciones más 
utilizadas para representarlas son la ecuación de 
van Genuchten (1980), y la ecuación de Brooks 
y Corey (1964), respectivamente: 


Ed a] (2) 
. y 
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0-0, Y 
K0)ox, [Pe (8) 


donde 6, es el contenido de humedad a satu- 
ración, generalmente asimilada a la porosidad 
total del suelo (p); 0, el contenido de humedad 
residual; Y, la presión; Y y Un valor característi- 
co de la presión; K, la conductividad hidráulica 
saturada; n, m y n son parámetros de forma 
adimensionales y positivos. 

Los parámetros m y n son relacionados 
a partir de la restricción propuesta por van 
Genuchten (1980) al modelo de la conductivi- 
dad hidráulica de Burdine (1953): m = 1-2/n. 
El parámetro n es relacionado con 1 y n, y la 
porosidad total mediante (Fuentes, Zavala, éz 
Saucedo, 2009): 


n=2s ( + 1 (4) 
mn 

donde s = Dj! E, en la cual D, es la dimensión 
fractal del medio poroso y E =3 es la dimensión 
de Euclides del espacio físico, ligado a la poro- 
sidad total con la ecuación implícita siguiente: 


(1-4) +9” =1 (5) 


Existen algunos trabajos en la literatura que 
estudian la ecuación de Richards en dos dimen- 
siones: Zavala-Trejo, Fuentes-Ruíz y Saucedo- 
Rojas (2005) estudiaron el drenaje agrícola sub- 
terráneo, para lo cual utilizaron la ecuación de 
Richards bidimensional (2D). Saucedo, Zavala, 
Fuentes y Castanedo (2013) acoplaron la ecua- 
ción de Barré de Saint Venant unidimensional 
con la ecuación de Richards 2D para modelar el 
riego por melgas. López-Avendaño et al. (1997) 
analizaron la infiltración de agua en un dominio 
tipo surco utilizando la ecuación de Richards 
2D. Sin embargo, la solución que aplican para 
resolverlo es con la discretización de la derivada 
de manera convencional: agregando un factor 
de peso o parámetros de interpolación para el 
tiempo y el espacio. En este estudio se resuelve 
la ecuación de Richards 2D —aplicado al riego 
por surcos— usando dos esquemas: un esquema 
implícito para todo el perfil y un esquema mixto: 
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explícito en la zona no saturada e implícito en 
la saturada. El esquema implícito se resolvió 
con el gradiente conjugado para problemas no 
lineales usando en la superficie una función 
sinusoidal. El modelo propuesto está adaptado 
para simular la variación del tirante en el tiem- 
po (condición de Dirichlet) y en la zona seca 
del surco se puede imponer una condición de 
tipo Neumann, por lo que pueden aplicarse 
condiciones de evaporación o una función de 
transpiración de acuerdo con las necesidades 
que se requieran. 


Materiales y métodos 
Creación del dominio espacial 


El dominio espacial para resolver la ecuación 
de Richards se creó con las siguientes etapas: 
1) selección y discretización de las fronteras del 
dominio; 2) relleno de las regiones interiores con 
puntos aleatorios distribuidos uniformemente, 
y 3) triangulación con el algoritmo de Delaunay 
(1934). Los dos tipos de fronteras existentes son 
exteriores (donde los puntos de interés están 
en el interior) e interiores (donde los puntos de 
interés están en el exterior). 

Para modelar la infiltración en los surcos se 
supuso una frontera cuadrada y el lado superior 
se sustituyó por una función periódica, en este 
caso, se optó por una función sinusoidal, 
Ax) = 15 sin (nx /45), que da un surco de 30 cm 


de altura y 90 cm de lomo a lomo (figura 1). Una 
vez que se obtuvo la frontera, el siguiente paso 
fue obtener la discretización de la misma. 

Para rellenar el interior primero se obtuvo 
el cuadrado más pequeño que contiene a todos 
los puntos frontera. Una vez obtenido este 
cuadrado C: [a, b] x [c, d] se buscó el punto P, 
usando una distribución uniforme dentro del 
cuadrado C. Para saber si P, estaba dentro de 
la frontera, se encontró la suma de los ángulos 
desde P, de cada par de puntos contiguos de la 
frontera (figuras 2 y 3). 

Para un punto en el interior de una frontera, 
la suma de los ángulos mencionados es de +360 
grados y para un punto en el exterior de la 
frontera la suma es de cero grados. Al descartar 
los puntos exteriores queda un dominio (figura 
3(a)). La gran cantidad de puntos en la figura 
3(a) que se generaron en este primer paso es 
porque se requiere que las regiones interiores 
estén lo más densas posibles. Posteriormente, 
se procedió a eliminar puntos que estuvieran a 
una distancia que se consideró conveniente, por 
ejemplo, todos aquellos que estuvieron a una 
distancia menor de 5 mm uno del otro (figura 
3(b)), para finalmente proceder a la unión de 
los puntos restantes, obteniéndose como resul- 
tado una triangulación del dominio, como el 
mostrado en la figura 3(c). El dominio con más 
puntos tendrá una mayor cantidad de nodos, lo 
que implica un costo computacional mayor para 
encontrar la solución. 


Figura 1. Izquierda: frontera con la parte superior parecida a un surco de riego. La parte superior f es una función sinusoidal 


trasladada y escalada para que concuerde con la parte superior de los lados a y c. Derecha: discretización de la frontera 


izquierda, puntos frontera. 
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. . . . Puntos frontera 


+ "Ángulo cubierto a cada par de 
puntos frontera contiguos 


Figura 2. Punto Pk dentro de una frontera tipo “surco”. Obsérvese que para un punto que está en el interior de una frontera, el 
ángulo cubierto es de 360 grados. 


Figura 3 (a): dominio tipo “surco” con 10,000 puntos interiores; (b): dominio (1) después de eliminar los puntos que tienen una 


distancia menor a 5 mm; (c): triangulación del dominio (2). 


No obstante que la función que se utilizó que HO _2 xp 00) 
para representar un surco es periódica y que dt  0x 9x 
con una mitad sería suficiente para modelar a oy (x) ok [w (x)] 
la infiltración, en este trabajo se hace toda la ea K [w(k)] E e (6) 
ejemplificación usando tres surcos; esto debido 
a que se realiza una simulación con cargas dife- 
rentes en cada uno de los surcos con la finalidad En esta ecuación, “(k) es una función de tres 
de observar el comportamiento del perfil de variables W(x, z, t). La ecuación (1) con derivada 


bumedecioiente: temporal discreta adquiere las siguientes dos 


formas (además de otras, dependiendo del 


esquema y precisión elegidos): 
Solución en diferencias finitas para la 


ecuación de Richards 2D 


Yw -W, 
c(u,) 2-2 1x(w)% 
t, > tos di dx 
La ecuación (1) evaluada en algún punto t en el 
. y La y 0 vw, 9K (w, ) 
tiempo se escribe de la siguiente forma: +=|K(w,) EA (7) 
0Z Oz 0Z 
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les we, 
el y tt (Ku) Ox 
re data 6) 
OZ OZ OZ 


La ecuación (7) corresponde al modelo ex- 
plícito y la ecuación (8) corresponde al modelo 
implícito. Así, se define el lado derecho de la 
ecuación (8) como el operador Op(Y ); las ecua- 
ciones (7) y (8) toman una forma más simple: 


c(w,) Lee - op (w,) (9) 
Eos =t, 

cu) Pra = Op (W,) (10) 
t, tu 


Cuando C[W(X)] +0, la solución que se obtie- 
ne es una solución recursiva: 


> a op (9) + (11) 


c(w) 
Sin embargo, en la zona saturada C(W,) =0, 
por lo que en esos puntos no puede aplicarse la 
ecuación (11). Para estos casos hay que optimi- 
zar la ecuación (10). El algoritmo que se utilizó 
en este trabajo es el del gradiente conjugado 
para problemas no lineales, también llamado 
“gradiente conjugado no lineal”. El método 
requiere que se defina una función a minimi- 
zar, en nuestro caso requerimos que se cumpla 
la ecuación (10), por lo que naturalmente una 
función de error es la siguiente: 
E(W,)= elo 


ya (12) 


1 — 
La E Op LY.) 


El gradiente conjugado no lineal para 
Richards 2D 


Siguiendo el procedimiento descrito por Dai y 
Yuan (1999), a partir de la ecuación (12), defini- 
da como la función de error, y con el operador 
definido a partir de la ecuación (8), el procedi- 
miento consta de los siguientes pasos: 1) se 


escogió x, = W(k); 2) Ax, = -V_E(x,); 3) a, = arg 
min [E(x, + a. : Ax,)]; 4) x, =x, +0, * Ax, y 5) s, = 
Ax,. Del paso 1 al paso 5 es la parte donde inicia 
X los siguientes pasos (6 - 10) son iterativos y el 
algoritmo termina cuando E(x,) es pequeño: 6) 
Ax, =-V,E(x,);7) B=(Ax7,: Ax,)/(Ax7,, - Ax, ,); 8) 
s,= Ax, +B>s,,;9 a, =arg min [E(x, +0: Ax )] 
y 10)x,,,=x,+a, *s . Por último se escoge x 
como la solución a las presiones en el tiempo 


n+1 n+1 


t. Aunque es posible aplicar el algoritmo para 
todo el vector de presiones, por razones de 
rendimiento los cambios se realizaron sólo en 
la zona saturada, debido a que, en la zona no sa- 
turada, es posible aplicar la regla de recurrencia. 
La condición inicial para x, es W(k) en la zona 
no saturada con el método explícito y W(k-1) en 
la zona saturada. 


Método para resolver el operador Op[W(k)] 


Para simplificar la notación se hicieron las 
siguientes simplificaciones: se define Y, como 
la presión en el nodo í en cualquier tiempo; 
la coordenada espacial del nodo ¡ se define 
como C(x,z,), donde x, es la coordenada en la 
longitud y z, en la profundidad; C, es la coor- 
denada espacial del j-ésimo nodo relacionado 
con el nodo i, es decir Co = Es = (X, ¿,2, m, es 
la cantidad de nodos relacionados con el nodo 
i, y Pi es la terna (x,z,W.) = (C,W.). Utilizando 
derivación de orden uno, la derivada en C, es 
estimada derivando el plano a:x+b-:z+c, 
que mejor ajusta a los puntos (P,, Puy Esp K, 
P,). Una forma de aproximar dicho plano es 
utilizar la matriz pseudoinversa definida por C,, 
y resolviendo para a, b y c; los puntos (P,, Pus 
Pape 
ya:x +b 7 +c=W ya-:x, +b:z2 +c=W 


) satisfacen a - x + b + z +c, por lo que 


mi 
ij 
Las ecuaciones anteriores definen el siguiente 
sistema de ecuaciones: 


x Z yw 

i i 1 1 

a Xx. Z, 1 a yw 
Albl=| a Alpe] 2 | (13) 

C ; : R C : 

2% z al y 
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La matriz A de la ecuación (13) no necesaria- 
mente es invertible; sin embargo, partiendo de 
la hipótesis que todo sistema se puede optimizar 
(Penrose, 1955), queda definida como: 


(14) 


NA 
ll 
> 
+ 
= 


donde A* es la matriz pseudoinversa de A. 
Con los valores a, b y c obtenidos, es posible 
derivar hacia cualquier dirección, en particular, 
hacia x y z. De manera similar a la aproximación 
lineal, utilizando derivación de orden dos, la 
derivada en C, se estimó derivando la superficie 


a: 2+b:x"y+c:y+d:x+e:y+fque mejor 
). Una 


mi 


ajustara a los puntos (P,, E E K,P 
forma de aproximar dicha superficie es utilizar 
la matriz pseudoinversa definida por C,, y re- 
solviendo para a, b, c, d, e, y f; así se obtuvieron 
los puntos (P, Pa Paj K, Pal que satisfacen 
a:x+b:x-y+c:y+d:x+e-y+f porlo que 
también satisfacen 4: 9+b+x,: y, +Cc:yi+d: 
x +e y +Hf=W ya: xi +b-x, y, +0 ys +d 


«x, +e:y, +f=W,. Estas ecuaciones quedan 
ij i,] ij 
definidas de la siguiente manera: 


Xx; Z; 1 w, 
a ss z a 
ajo 1 1, || “2 [(15) 
C ' / C 
x z 1 


im; im im; 


donde aplicando la misma consideración que 
en la solución de la ecuación (14), la matriz A 
queda definida como: 


q 
b w, 
y 
C + il 
=A"- ' (16) 
d : 
e Mo 


donde A* es la matriz pseudoinversa de A. 
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Al final de este procedimiento, con los 
valores a, b, c, d, e y f obtenidos, se calculó la 
derivada. 


Cálculo de la lámina infiltrada 


El volumen infiltrado por unidad de longitud 
de cauce en la unidad de tiempo (4) se calculó 
con la siguiente expresión (Fuentes, De León- 
Mojarro, € Hernández-Saucedo, 2012): 
o JA 
eee A(xt)= $9,(x,y,2,4)dP, (17) 
P 


donde q, (x,y,z,t) es el caudal de infiltración por 
unidad de superficie de canal o caudal unitario 
y P, es el perímetro mojado en cada surco. 


Aplicación 


Comparación entre el esquema implícito y el 
esquema mixto 


Los dos esquemas (implícito y mixto) se com- 
pararon. Los valores que se utilizaron fueron 
del suelo Montecillo (Fuentes, 1992; Saucedo, 
Zavala, é Fuentes, 2011; Saucedo et al., 2013), los 
cuales son 6, = 0.0000 cm? cm”, O, = 0.4865 cm? 
cm”?, K_ = 1.8400 cm h”, Y, = -45.9000 cm, 
m = 0.1258 y una profundidad L = 90 cm. En esta 
simulación se utilizó una distancia promedio 
entre nodos en ambos esquemas de Az=1 cm, y 
una discretización temporal de At =0.01 h para 
el método mixto, mientras que para el implícito 
se usó un At =0.00128 h. Los resultados para 240 
min (figura 4) indican que no hay diferencias 
entre ambos esquemas, pero el esquema mixto 
requiere una discretización temporal mucho 
más fina que el esquema implícito. Sin embargo, 
el esquema implícito requiere más cálculos para 
cada tiempo, lo que hace la simulación de un 
evento de infiltración un poco más lenta. 

El tiempo de cómputo requerido para si- 
mular los 240 min, usando una computadora 
de 6 GB en RAM, para el método mixto fue de 
18 min, mientras que para el método implícito 
fue de 24 min aproximadamente, es decir, el 
método implícito requirió de un 30% de tiempo 
adicional que el método mixto, sin embargo, con 
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3 3 
o TE o A 
me cm? 


180 cm 


90 cm 


Figura 4. Comparación a los 240 minutos entre el esquema implícito (líneas blancas) y el esquema mixto (líneas negras). 


los recursos computacionales de hoy en día, el 
tiempo de simulación puede acortarse un poco 
más. 

Es importante mencionar que este tiempo 
de cómputo corresponde a la discretización de 
los tres surcos; si esta simulación se realizara 
únicamente en la mitad de un surco, los tiempos 
para obtener la solución se reducirían con la 
misma discretización espacial y temporal a 4.5 
y 6 minutos, respectivamente. El resultado de 
las curvas de infiltración de las soluciones son 
prácticamente las mismas. 


Simulación de la infiltración en surcos en 2D: 
tirantes constantes en los surcos 


Tomando como base las características hidro- 
dinámicas del suelo Montecillo se realizó una 
simulación del proceso de infiltración en tres 
surcos, en donde se aplicó un tirante de 10 cm 
en la primera y la tercera hora de infiltración, y 
en la segunda y cuarta hora se aplicaron 5 cm 
de tirante. Se trazaron curvas con los siguientes 
contenidos volumétricos de agua: 0.47 cm' cm”, 
0.43 cm? cm?, 0.39 cm? cm?, 0.35 cm? cm”, 0.31 
cm? cm”. En la figura 5 se muestran dos perfiles 
de humedad correspondientes a la primera hora 
de infiltración, en el minuto 60 el tirante cambia 
de 10 a 5 cm. Como puede verse, el frente de 
humedecimiento avanza más rápido en la pri- 
mera hora (figura 5) que en la segunda (figura 
6) debido a que en la segunda hora se tiene un 


tirante de 5 cm desde la base del surco, por lo 
que, en dicha hora, los efectos de la capilaridad 
aportan más que los efectos de presión. El efecto 
inverso puede observarse en la imagen superior 
de la figura 7, en el cual la presión hacia abajo es 
mayor que en los minutos anteriores, debido a 
que el tirante cambió a 10 cm en el minuto 120. 

La figura 6 muestra dos perfiles de hume- 
dad correspondientes a la segunda hora de 
infiltración y la figura 7 muestra dos perfiles de 
humedad correspondientes a la tercera y cuarta 
hora. Como puede verse en la serie de imágenes, 
cuando se hace el cambio de tirante de 10 a5 cm 
en el surco se favorece el humedecimiento en 
los lomos del surco, mientras que con un tirante 
mayor la redistribución del perfil de humedad 
es mayor en la base del surco que en el lomo. 
Con un tiempo de simulación de 60 minutos se 
puede apreciar que ya se ha logrado llevar el 
contenido de humedad a 0.31 cm? cm? en una 
profundidad de 35 cm y los lomos del surco. 

El contenido de humedad en el surco a los 
180 minutos de simulación puede verse que está 
prácticamente saturado, es decir, tiene una can- 
tidad de agua ligeramente mayor a 0.47 cm? cm? 
y la máxima cantidad de agua que puede tener 
es de 0.4865 cm? cm”. Después de los 135 min, 
la parte superior de los surcos está saturada, y 
a partir de ese momento, para cualquier conte- 
nido de humedad en el perfil del suelo tiende a 
ser horizontal, que se asemeja a una infiltración 
constante (figura 7). 
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Añ 


30 cm 


4 


90 cm 


¿ cmó cmd 
== 035 e. Eos. 


cm 


180 cm 


Figura 5. Primera hora de infiltración con un tirante de 10 cm desde la base del surco. A los 60 minutos el tirante cambia a 5 cm. 
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Figura 6. Distribución de la humedad durante la segunda y hora. El tirante de 5 cm es constante. 
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Figura 7. Distribución de la humedad en la tercera y cuarta hora. 


Simulación de la infiltración en surcos en 2D: 
tirantes diferentes en los surcos 


Usando las mismas características hidrodinámi- 
cas del suelo Montecillo del inciso anterior se 
realizó una simulación en tres surcos, en donde 
en cada uno de ellos las cargas ahora fueron 
diferentes: 3, 9 y 6 cm, de izquierda a derecha. 
Los resultados de la simulación mostraron que 
la redistribución de humedad en los lomos del 
surco y en el perfil del suelo durante el proceso 
de infiltración es mejor en el área entre los sur- 
cos con cargas 9 y 6 cm (área 1), que el que está 
entre los surcos con 3 y 9 cm de carga (área 2) 
(figuras 8 y 9). 

La distribución de la humedad en los lomos 
de los surcos, la unión de los bulbos de hume- 
decimiento y área de mojado es mejor en la 
primera área que en la segunda. En el área 1, el 
surco que tiene 9 cm de carga tiene que aportar 
mayor cantidad de agua al surco de la izquierda 
(carga de 3 cm), por lo que la distribución de la 


humedad hacia capas inferiores es menor que 
entre los surcos de 9 y 6 cm. 

Además, el contenido de humedad en el 
lomo del surco con carga de 9 cm es mayor que 
en el lomo del surco con 3 cm; en este último, el 
principal problema es que se retrasa el mojado 
hacia las capas inferiores (figura 9). Este pro- 
blema es muy común en los eventos de riego 
por gravedad, donde el caudal aplicado en cada 
surco no es homogéneo debido principalmente 
a que los sifones no están calibrados, que hay 
fallas en la apertura de compuertas o, como en 
la mayoría de los casos, que el gasto propor- 
cionado en cada surco es a criterio del regador. 
Lo anterior trae consigo varias complicaciones, 
entre la más importante está que la lámina 
aplicada en cada surco no es homogénea, lo que 
conduce a que en algunas partes de la parcela se 
aplique más agua de la que necesita o, en el peor 
de los casos, que se aplique una cantidad menor. 

Por tratarse un problema en dos dimensio- 
nes, la lámina infiltrada es ahora representada 
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Figura 8. Distribución de la humedad en la primera y segunda hora. 
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Figura 9. Distribución de la humedad en la tercera y cuarta hora. 
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como el área infiltrada, que si se multiplica por 
una longitud se obtiene el volumen infiltrado. 
Así, el área infiltrada (cm?) se calculó en cada 
surco con la ecuación (17) y al hacer las com- 
paraciones (figura 10) se puede ver que para el 
tiempo de 240 min, los surcos de 9, 6 y 3 cm de 
tirante infiltran en promedio 512, 467 y 401 cm?, 
respectivamente. Es decir, se estarían infiltrando 
110 cm? menos en el surco con carga de 3 cm; 
además de que en los surcos con mayor tirante 
tenderán, por diferencia en los tirantes, a contri- 
buir al humedecimiento de los surcos con menor 
tirante, lo que también trae consigo una dismi- 
nución en el patrón de mojado en cada uno de 
los surcos que les corresponde humedecer. 


Conclusiones 


El fenómeno de infiltración en un dominio bidi- 
mensional se modeló y resolvió mediante el mé- 
todo del gradiente conjugado. La naturaleza de 
la ecuación de Richards hizo que no fuera sufi- 
ciente el método explícito —a diferencia de otras 
ecuaciones no lineales—, por lo que se combinó 


con el método implícito. La metodología pre- 
sentada mostró buenos resultados, además de 
que muestra la bondad de incorporar un tirante 
variable en el tiempo. 

La arbitrariedad en la posición de los 
puntos del dominio espacial nos permitió 
crear cualquier discretización espacial; sin 
embargo, esto ocasiona que las derivadas 
tengan que ser resueltas usando la superficie 
que más se acerque a una región de puntos. 
Con el uso de esta metodología, en la zona 
seca del surco se puede imponer una condición 
de tipo Neumann, por lo que pueden aplicarse 
condiciones de evaporación o una función de 
transpiración de acuerdo con las necesidades 
que se requieren. 

Finalmente, la solución utilizada puede ayu- 
dar a tomar decisiones sobre el tiempo de riego 
y la profundidad a la cual se desea mojar, ya 
que en las simulaciones realizadas con tiempos 
máximos de 4 h se pudo ver que la profundidad 
saturada alcanzó los 50 cm, tiempo suficiente 
como para mojar la zona radicular de la mayoría 
de los cultivos. Sin embargo, en la práctica, 
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Figura 10. Evolución de la lámina infiltrada en surcos con carga variable. 
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estos tiempos son mucho mayores, lo que trae 
consigo la aplicación de láminas excesivas y 
bajas eficiencias de aplicación. 
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